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Abstract: Non-erodible elements such as stones and vegetation are key to controlling wind erosion and 
dust emission in drylands. Stony deserts are widely distributed in the Gobi Desert, but the effect of 
stones on wind erosion and dust emission have not been well studied, except under artificial conditions. 
In this study, we evaluated the effect of stones on wind erosion and dust emission by measuring the 
sand saltation threshold in a stony desert in Tsogt-Ovoo in the Gobi Desert, Mongolia, under natural 
surface conditions during sand and dust storms. We quantified the amount of stones by measuring the 
roughness density, and determined the threshold friction velocity for sand saltation by measuring wind 
speed and sand saltation count. Our results showed that the threshold friction velocity increased with 
the roughness density of stones. In the northern part of the study area, where neither a surface crust 
nor vegetation was observed, the roughness density of stones was 0.000 in a topographic depression 
(TD), 0.050 on a northern slope (N.SL), and 0.160 on the northern mountain (N.MT). The mean 
threshold friction velocity values were 0.23, 0.41, and 0.57 m/s at the TD, N.SL, and N.MT sites, 
respectively. In the southern part of the study area, the roughness density values of stones were 0.000 
and 0.070-0.320 at the TD and southern slope sites, respectively, and the mean threshold friction 
velocities were 0.23 and 0.45—-0.71 m/s, respectively. We further compared the observed threshold 
friction velocities with simulated threshold friction velocities using Raupach's theoretical roughness 
correction and the measured roughness density values, and found that Raupach's roughness correction 
worked very well in the simulation of threshold friction velocity in the stony desert. This means that the 
results of our stone measurement can be applied to a numerical dust model. 
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1 Introduction 


Sand and dust storms can cause human and animal mortality (Mu et al., 2013) and health 
deterioration (Onishi et al., 2015), as well as agricultural production losses (Yang et al., 2002). 
Additionally, aeolian dust particles that originate from such storms have an impact on climate, 
though the scale of that impact is still under debate (e.g., Boucher et al., 2013). To avoid 
damage from storms and to evaluate their effects, a high-precision three-dimensional numerical 
dust cycle model is needed because of the large uncertainty in the amount of dust emission of 
current models, which shows an 11.1-fold difference in global dust emission (Wu et al., 2020). 
Many factors influence erodibility and dust emission. Factors affecting the susceptibility of soil 
and land surfaces to wind erosion include the presence and coverage of stones, soil particle size, 
the presence of a soil crust, soil wetness, and the amounts of snow cover and vegetation 
(Kurosaki et al., 2011a). The effects of these factors on the occurrence of sand and dust storms 
in East Asia, one of the world's major dust sources, have been evaluated by using an erodibility 
index, such as the threshold wind friction velocity, and dust occurrence ratios (Kurosaki et al., 
2011a, b; Wu et al., 2016, 2020). However, the effect of stones has not been evaluated under 
natural conditions in East Asia. The Gobi Desert in northern China and southern Mongolia, 
along with the Taklimakan Desert in Northwest China, is a major dust source in East Asia (e.g., 
Natsagdorj et al., 2003; Wang et al., 2008; Chen et al., 2017). The Gobi Desert is characterized 
by widely distributed stone pavements (Cooke, 1970). An area covered by a stone pavement is 
referred to as a "gobi", which means "stony desert" in Chinese, but in Mongolian, "gobi" is used 
for a "desert with sparse vegetation, barren soil, and scarce surface water". 

Aeolian sand saltation is a key dust emission process (Greeley and Iversen, 1985; Shao et al., 
1993). The threshold value of wind friction velocity for sand saltation depends on many 
parameters, including soil particle size (Gillette et al., 1980; Shao and Lu, 2000; Edwards and 
Namikas, 2015), the presence of non-erodible elements (Gillette, 1978; Marticorena et al., 2006; 
Tan et al., 2019), soil moisture (Ishizuka et al., 2009), soil crust (Gillette et al., 1982; Marticorena 
et al., 1997a), and other soil properties (Gillette, 1988). Non-erodible roughness elements, which 
include vegetation, stones, and soil crusts, influence the major physical processes involved in 
wind erosion (Callot et al., 2000) through two effects: the coverage effect and the shear stress 
(i.e., drag) partitioning effect (Marshall, 1971; Iversen et al., 1990). In quantitative terms, 
increased coverage of non-erodible elements, by reducing the total erodible area, reduces the 
sand saltation amount. Similarly, the partitioning of the total shear stress due to wind into the 
force exerted on the roughness elements and that exerted on intervening surfaces (Schlichting, 
1936) reduces the wind shear stress that reaches the soil surface (Gillies et al., 2007) and 
increases the sand saltation threshold. Li et al. (2020) have shown that variation of sand saltation 
threshold is due to the different shear stresses caused by the presence of non-erodible elements. 

Several theoretical models have been developed to describe the partitioning of shear stress 
between the surface and roughness elements (Raupach, 1992; Raupach et al., 1993; Marticorena 
and Bergametti, 1995). Specifically, these models describe the effect of roughness elements on 
the threshold friction velocity, which is the minimum friction velocity required to initiate 
movement of sand particles. In dust models, the frequency of sand and dust events is determined 
by how often the wind velocity exceeds the threshold friction velocity, and their intensity 
depends on by how much the threshold is exceeded (Callot et al., 2000). The model developed 
by Raupach et al. (1993) accounts for the effects of roughness elements and their coverage on 
the threshold friction velocity (Webb et al., 2020) by using the roughness density, which is 
defined here as the sum of the element frontal areas (i.e., facing the wind direction) divided by 
the ground surface area, and constant parameters such as the basal-to-frontal area ratio of non- 
erodible elements. In general, the model by Raupach et al. (1993) performs well with respect to 
the results of wind tunnel experiments (Gillies et al., 2000), field experiments conducted with 
artificial cylindrical roughness elements (Gillies et al., 2007), and field experiments conducted 
in vegetated areas (e.g., King et al., 2005; Gillette et al., 2006). 

To investigate the effects of two types of roughness elements, i.e., solid elements (e.g., 
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gravel, stone) and flexible elements (e.g., vegetation), on the threshold friction velocity, 
Marticorena et al. (1997b) developed a double shear stress partition model based on the single 
shear stress partition model of Marticorena et al. (1997a), and Darmenova et al. (2009) 
developed a double shear stress partition model based on the model of Raupach et al. (1993). 
Darmenova et al. (2009) used their model to simulate dust transport over East Asia, as did Xi 
and Sokolik (2015), Foroutan et al. (2017), and Lee and Park (2020) using the same model. All 
of the above-mentioned model studies employed representative roughness density of stones for 
each land cover type because data on the roughness density of a natural stone surface covering 
wide areas of the world were not available. For example, the minimum roughness density of 
solid elements is constant in barren and sparsely vegetated land, which is the dominant land 
cover type in our study area (Tsogt-Ovoo in the Gobi Desert, Mongolia). We question whether 
accurate results can be obtained by using a constant representative roughness density for natural 
stones to estimate the threshold friction velocity in this study. 

The effects of stones on wind erosion in the Gobi Desert of northern China have been studied 
mainly by field wind tunnel measurements (e.g., Tan et al., 2013; Zhang et al., 2014; Zhang et 
al., 2015; Zhang et al., 2016). In these above-mentioned studies, stone coverage was strongly 
correlated with the threshold friction velocity, and the sand flux decreased exponentially with 
stone coverage. However, these studies did not consider the effect of the roughness density of 
stones on the threshold velocity. Recently, Tan et al. (2019) conducted a wind tunnel experiment 
in the Gobi Desert for aeolian saltation to determine the threshold wind velocities for artificial 
stony surfaces with different roughness densities and revealed that the threshold friction 
velocity for sand saltation increased with the roughness density of stones. Their results agree 
well with the model of Raupach et al. (1993). 

Field-based, wind tunnel, and model simulation studies have all shown that the aerodynamic 
roughness length and roughness density of non-erodible elements are important parameters for 
assessing the effect of surface roughness elements on the threshold friction velocity in wind 
erosion models. However, roughness elements are generally described only by their roughness 
density (Marticorena et al., 1997b). There is a good relation between the ratio of aerodynamic 
roughness length to the mean height of roughness elements and roughness density (Marshall, 
1971; Lancaster and Baas, 1998; King et al., 2005; Marticorena et al., 2006; Hébrard et al., 
2012). Thus, for natural stones, if the roughness density and the mean height are well quantified, 
then the aerodynamic roughness length can be obtained by using the aforementioned relation 
and could be used in wind erosion models. Therefore, in this study, we measured the roughness 
density of stones and the observed sand saltation thresholds in Tsogt-Ovoo in the Gobi Desert, 
Mongolia. We also used the model of Raupach et al. (1993), which includes a roughness 
correction, to simulate threshold friction velocity using the measured roughness density of 
stones. The main objective of this research was to evaluate the effects of stones on the threshold 
friction velocity for sand saltation. The specific objectives were as follows: (1) to quantify 
roughness density for different coverages, numbers, and sizes of stones; (2) to calculate the 
threshold friction velocity using observed sand saltation count data and meteorological data; and 
(3) to validate the theoretical model with roughness correction by simulating threshold friction 
velocity using the measured roughness density of stones. 


2 Study area and methodology 
2.1 Study area 


Sand and dust storms are a common phenomenon in Asian drylands such as the Gobi Desert and 
Taklimakan Desert. Our study sites are located in Tsogt-Ovoo Village (44°23'04"N, 
105°16'59”E), South Gobi Province, Mongolia, which is in the northern part of the Gobi Desert 
(Ishizuka et al., 2012). An analysis of meteorological observatory data has shown that Tsogt- 
Ovoo experiences the highest frequency of dust storms in East Asia (Kurosaki and Mikami, 
2007). The area around Tsogt-Ovoo includes a topographic depression that extends 
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approximately 10 km in the northeast-southwest direction and 20 km in the northwest-southeast 
direction (Fig. 1). Further, the vegetation coverage in Tsogt-Ovoo differs depending on the 
topography (Gantsetseg et al., 2017). The mean annual precipitation during 1975-2019 was 
100.4 mm. More than two-thirds of the annual precipitation occurs from June to August. The 
annual mean temperature was 3.8°C in 2011 with a maximum temperature of 37.1°C in August 
and a minimum temperature of —31.3°C in January (Ishizuka et al., 2012). 
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Fig. 1 Topographic map showing the locations of the observation sites around Tsogt-Ovoo. Stars indicate sites 
where observations were conducted year-round: The Main site (labeled 1; see Table 1) was established in 
March 2012, and sites Sub14A and Sub14B (labeled 2 and 3, respectively) were established in March 2014. 
Temporary observation sites are indicated by triangles. Observations were conducted during 1-5 May 2018 at 
sites Sub18A—Sub18D (labeled 4-7, respectively) and from 27 April to 8 May 2019 at sites Subl19A—Sub19D 
(labeled 8-11, respectively). Wind speed, wind direction, and sand saltation count were observed at all sites, 
whereas soil moisture was observed only at the year-round sites. 


The sites are located in five different topographic positions, designated southern mountain 
(S.MT; Sub18A site), southern slope (S.SL; Sub14B, Sub19A, and Sub19B sites), topographic 
depression (TD; Sub14A, Sub18D, Sub19C, and Sub19D sites), northern slope (N.SL; Main and 
Sub18B sites), and northern mountain (N.MT; Sub18C site) (Table 1). It should be noted that 
the two digits in the site name indicate the installation year. We chose the temporary site 
locations on the basis of land surface features, especially the size and amount of stones. 

Figure 2 shows land surface conditions in front of the sand saltation count sensor at the 
different sites. There was no green or dead vegetation at the Main, Sub14A, and Sub14B sites in 
spring 2018, but green or dead vegetation was observed at those sites in spring 2019. There was 
a hard crust at the Sub14A, Sub18A, and Sub18D sites in both years. There was a soft crust on 
sandy soil at the Sub19C site. Further, at the site Sub19D, there was a sand sheet with no 
vegetation, crust, or stones. 


2.2 Experimental design 


Because land surface characteristics differed from site to site in our study area, we compared the 
simulated threshold friction velocity at each study site with the observed value to validate the 
theoretical model. 


202108.00011v1 


chinaXiv 


ChinaXivA (ERAT! 


Batjargal BUYANTOGTOKH et al.: Effect of stones on the sand saltation threshold during... 


Table 1 Site locations and type, measurement periods, topographic positions, and surface characteristics 


No. Site Geographical coordinates Lype of Obar anon Topographic Surtace, 
name site period position characteristic 
1-5 May 2018 Stone 
1 Main  44°23'04"N  105°16'59"E Year-round 27 April-8 May N.SL Stone and 
2019 vegetation 
1-5 May 2018 Hard crust 
2  Subl4A  44°23'37"N  105°13'24"E Year-round 27 April-8 May TD Gerad 
2019 vegetation 
1-5 May 2018 Stone 
3  Subl4B 44°23'32"N  105°13'12"E Year-round 27 April-8 May S.SL Stoneand 
2019 vegetation 
4 Sub18A 44°22'23"N 105°11'28"E Temporary 1-5 May 2018 S.MT Stone 
5 Sub18B 44°24'26'"N 105°14'01"E Temporary 1-5 May 2018 N.SL Stone 
6 Sub18C 44°27'32'N 105°16'48"E Temporary 1-5 May 2018 N.MT Stone 
7 Sub18D 44°21'47"N 105°14'S51"E Temporary 1-5 May 2018 TD Hard crust 
8 Subl9A  44°21'52"N  105°14'04"E Temporary 7” E May S.SL Stone 
9  Subl9B  44°22'02"N  105°14'27"E Temporary 7” A May S.SL Stone 
10 Subl9C  44°22'29"N  105°14'54"E Temporary 77 APril-8 May TD Sand and Sot 
2019 crust 
11 Subl9D  44°22'39"N  105°14'59"E Temporary 77 APril-8 May TD Sand 


2019 
Note: S.MT, southern mountain; S.SL, southern slope; TD, topographic depression; N.SL, northern slope; N.MT, northern mountain. 
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Fig. 2 Land surface conditions in front of the saltation count sensor at 14 observation sites. Each site name is 
shown at the bottom of each photo with the topographic position (MT, mountain; SL, slope; TD, topographic 
depression) and observation year in parentheses. The diameter of the mouth of the sand catcher was 3.1 cm, and 
the diameter of the sand saltation count sensor was 1.2 cm. 


The workflow for this study is shown in Figure 3. Specifically, image analysis was used to 
directly measure the height, width, and length of 105 stones at each site; these measurements 
were used together with the stone coverage to calculate the roughness density of stones. The 
obtained roughness density was used in a roughness correction parameterization to simulate the 
threshold friction velocity at each site and to determine the roughness length of the stony sites. 
The observed wind speed and roughness length were used to calculate the friction velocity. The 
obtained friction velocity and the observed sand saltation count (hereafter, saltation count) were 
used to determine the instantaneous threshold friction velocity. The instantaneous threshold 
values were used to calculate the mean threshold friction velocity during the observation period 
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shown in Table 1. Finally, the observed mean was used to validate the simulated threshold 
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Fig. 3 Workflow showing the observation and model simulation of the threshold friction velocity in this study. 
As is the roughness density of stones, and f(/s) is a correction functions for a roughness element. 


2.3 Observations of meteorological parameters, saltation count, and soil texture 


We measured the spatial coverage and size (i.e., height, width, and length) of stones, 
meteorological parameters such as wind speed and wind direction, saltation count, and soil 
texture at each site. Layout of quadrats where land surface photographs were taken for stone 
measurement and demonstration of anemometer and wind vane are shown in Figure 4. Wind 
speed, wind direction, and saltation count were measured at each site during 1-5 May 2018 and 
from 27 April to 8 May, 2019. Anemometers (S-WSA-M003, Onset Company, Bourne, MA, 
USA) and wind vanes (S-WDA-M003, Onset Company, Bourne, MA, USA) were installed at a 
height of 1.7 m at each site (Fig. 4b). Saltation count sensors (UD-101, Chuo Kosoku Ltd., 
Japan) (Udo, 2009) were installed at a height of 0.1 m at each site (Fig. 4b). Soil moisture 
sensors (S-SMC-M005, Onset Company, Bourne, MA, USA) were installed to 5 cm depth. The 
saltation count, wind speed, and wind direction were measured at 30 s intervals. The UD-101 
sensor was positioned facing the northwest (i.e., 315°), which is the prevailing direction of 
strong winds in the study area (Ishizuka et al., 2012). We used only the observational data 
obtained when the wind direction was in the northwest quarter, from 260° to 370° (or 10°). We 
collected soil samples from the surface layer (1-5 cm depth) at each site and analyzed the soil 
texture by the dispersion method (Faé et al., 2019). Soil moisture measurements were performed 
only at the Main, Sub14A, and Sub14B sites. 


2.4 Measurements of stones 


We used image analysis to obtain the breadth and height of stones and their number for 
calculating the roughness density of stones at each site. We took 21 photographs of the surface 
in 30 cmx50 cm quadrats (including the scale) at each site (Fig. 4a). We established three 21-m 
long observation lines extending northward (N), northwestward (NW), and westward (W) from 
a point in the southeast part of the site, where the measurement instruments for the saltation 
count and wind were installed. 
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D-101 sensor} 


Fig. 4 Layout of quadrats (gray rectangles) where land surface photographs were taken for stone measurement 
(a) and demonstration of anemometer and wind vane (b). The black dot in the left panel indicates the position 
of the UD-101 count sensor, which faced northwest. The UD-101 sensor was installed at the height of 0.1 m; 
the anemometer and wind vane were installed at the height of 1.7 m. 


2.4.1 Measurements of stone coverage, width, length, and number by image analysis 

We used ImageJ version 1.52p, developed by the U.S. National Institutes of Health to 
automatically count and measure the number and size of each stone, which are necessary to 
obtain the roughness density of stones (Mazzoli and Moriconi, 2014). Figure 5 shows the steps 
involved in image processing. 


Fig. 5 Demonstration of image processing steps. (a), cropped image with scale; (b), the color threshold for soil 
that is manually set; (c), image after application of the threshold for stones and the "Fill Holes" and "Watershed" 
functions to separate particles; (d), final image after application of the "Extended Particle Analyzer" function, 
which removes coarse sand particles. 


In the first step, the pixel size was measured by referring to the 30 cmx50 cm scale (see Fig. 
4a). After the scale was set, the image was cropped to exclude the scale (Fig. 5a). In the second 
step, the "Color Threshold" function was used to color the soil surface, while stones retained 
their original color. The color threshold was manually adjusted, but we detected some stones on 
the soil surface with a color similar to that of the soil by tracing their outlines by hand (Fig. 5b). 
In the third step, we converted the RGB image to an 8-bit grayscale image, setting the soil 
surface to 0 (black) and the stones to 255 (white). After this step, some gray areas were left. In 
the fourth step, these gray areas were divided into stone and soil automatically by using the 
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threshold method of Otsu (1979). However, some fine particles, too small to be stones, were 
left. Therefore, in the fifth step, the functions "Fill Holes" and "Watershed" were applied to 
separate these particles (Fig. 5c). Finally, in the sixth step, we revised colors of stones and 
background and removed fine particles with an area smaller than 0.12 cm’, corresponding to a 
diameter of 4.0 mm, by applying the "Extended Particle Analyzer" function (Fig. 5d). Then, we 
measured the width, length, and area of each stone and used the stone areas to estimate the stone 
coverage as the area fraction of stones (black areas in the final image) in each 30 cmx50 cm 
quadrat. The averaged stone coverage of the 21 images was considered to represent the stone 
coverage at the site. 


2.4.2 Roughness density of stones 


The roughness density, A; (1.e., lateral cover), is defined as the ratio of the total frontal area of all 
roughness elements to the total surface area, and the calculation formula is as follows: 


ns 2 bh, 
= — i=l 

ar a (1) 
where n is the number of roughness elements; b (cm) and A (cm) are the roughness element and 
height, respectively; s (cm?) is the mean frontal silhouette area of the obstacles present over the 
horizontal surface area S (cm?) (e.g., Marshall, 1971); and 7 is the stone. The frontal silhouette 
area of each stone was computed from the width and height of each stone. We determined the 
breadth of each stone by averaging the width and length of the stone obtained by the image 
analysis (Section 2.4.1). 

We estimated the height of each stone from the two-dimensional image by assuming that the 
breadth to height ratio at a given site was constant. To obtain the representative breadth to 
height ratio for each site, we measured the height, width, and length of 105 stones at the site in 
the field and used the averaged width and length values to calculate the breadth. We calculated 
the height of each stone in the two-dimensional images using the equation: h;=b,/r, where r is 
the representative breadth to height ratio. Finally, we used Equation 1 to calculate roughness 
density (As). The averaged roughness densities of the 21 images were considered to be the 
representative roughness density value at each site. 


2.5 Roughness length 


The aerodynamic roughness length is a function of the mean height, density, and shape of the 
solid elements. Several studies have empirically determined the relationships among 
aerodynamic roughness length, roughness mean height, and roughness density of natural desert 
surfaces, especially sparsely vegetated and stony surfaces (e.g., Lancaster and Baas, 1998; King 
et al., 2005; Marticorena et al., 2006; Hébrard et al., 2012). Foroutan et al. (2017) compiled 
aerodynamic roughness length, roughness density, and roughness mean height from data 
reported by several field and wind tunnel studies (King et al., 2005; Marticorena et al., 2006; 
Hébrard et al., 2012) and created a comprehensive database of aerodynamic roughness lengths. 
By fitting straight lines to the relationship between the ratio of aerodynamic roughness length to 
roughness mean height and roughness density, Foroutan et al. (2017) estimated the aerodynamic 
roughness length relevant to aeolian processes as follows: 


z, _ {0.960A,'"", 4 <0.2 
h |0.0834°, 4 20.2 


where zo (cm) is the aerodynamic roughness length; A (cm) is the roughness mean height; and A; 
is the roughness density. 

We used the mean stone height at each site and Equation 2 to determine the aerodynamic 
roughness length at the stony sites in this study. 

For relatively smooth surfaces, aerodynamic roughness length is approximately 1/30 of the 
mean particle diameter (Greeley and Iversen, 1985). The value of aerodynamic roughness length 
was used at the sandy sites or sites with crusts (Sub19C, Sub19D and Sub14A). 


(2) 
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2.6 Friction velocity and the threshold friction velocity estimated from observations 


2.6.1 Friction velocity 


The friction velocity, which represents the surface shear stress exerted by the wind, is a critical 
parameter in numerical dust emission schemes. It determines the onset and end of saltation, and 
it also affects the horizontal sand flux and the vertical dust flux. By assuming a logarithmic 
wind profile under neutral conditions, the friction velocity u» can be calculated as, 


= kU, J 


u, =——+_ , 3 
log(z/ zo) @) 
where k is the von Karman constant (k=0.41); zo (m/s) is the aerodynamic roughness length; z 
(m) is the observation height (z=1.7 m); and U7 (m/s) is the observed wind speed at the 
observation height of 1.7 m. Note that in Equation 3, we used the aerodynamic roughness length 


estimated from the roughness density obtained with Equation 2. 
2.6.2 Threshold friction velocity obtained from observations 


The threshold friction velocity for sand saltation is commonly defined as the minimum friction 
velocity required for the initiation of sand transport. Several methods have been employed to 
calculate this value (Barchyn and Hugenholtz, 2011), although generally the threshold wind 
velocity is actually calculated instead of the threshold friction velocity. We employed 
the'"Instantaneous Method" of Bagnold (1941), which identifies the friction velocity at which 
sand saltation begins and ends as the threshold friction velocity (Barchyn and Hugenholtz, 
2011). The instantaneous threshold friction velocity (w+: inst) was estimated every 30 s (Section 
2.3), and the mean threshold friction velocity (w+ mean) over the observation period (Table 1) was 
used for validating the simulated threshold friction velocities. Figure 6 presents an example of a 
friction velocity time series to illustrate the method used to determine the instantaneous 
threshold friction velocities. The black and blue gray lines show the friction velocity when sand 
saltation was observed or not observed, respectively. The red circles indicate the times at which 
sand saltation started and ended, and the friction velocities at those times were used as the 
instantaneous threshold friction velocities. The green line shows the friction velocity when the 
wind direction was between 10° and 260°. Because the sand saltation sensors (UD-101) faced 
northwest (i.e., 315°), the instantaneous friction velocities were not determined when the wind 
direction was between 10° and 260°, although the sand saltation was active. 
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Fig. 6 Time series of friction velocity during 28-29 May 2019 at the site Sub19D. Blue line, no sand transport; 
black line, active sand transport; w+: inst, instantaneous threshold friction velocity at the start and end of sand 
transport; w+: mean, mean threshold friction velocity (dashed horizontal line). Friction velocities obtained when 
the wind direction was between 10° and 260° (green line) were not used in the determination of the threshold 
friction velocity. Friction velocity was calculated using Equation 3 with aerodynamic roughness length of the 
threshold friction velocity. Friction velocity was calculated using Equation 3 with aerodynamic roughness 
length of 0.0012 m/s. Wind speed and wind direction (WD) were measured at 1.7 m height. UTC is the 
Coordinated Universal Time. 
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2.7 Parameterization of the threshold friction velocity in the numerical model 


The threshold friction velocity represents the minimum wind speed for the mobilization of soil 
grains. Under ideal conditions (i.e., no roughness elements, no soil moisture, no salt 
concentration, and no soil crust), the threshold friction velocity can be calculated by considering 
the balance among wind drag, gravity, and interparticle cohesion (Shao and Lu, 2000). The 
threshold friction velocity (u) for an arbitrary surface that is affected by roughness elements 
(shear stress partitioning), soil moisture, salt concentration, and soil crust is expressed as 
follows: 

Ue, =(d,A, W815) =U dS ASWI a) (4) 
where ds (mm) is the sand particle size; uo (m/s) is the threshold friction velocity of sand 
particles of size under ideal conditions; 2, w (m?/m°), sc, and ser are the roughness density of the 
roughness elements, soil moisture, salt concentration, and soil crust, respectively; and AA), fw), 
fsc), and fsc) are the correction functions for surface roughness elements, soil moisture, salt 
concentration, and soil crust, respectively (Shao et al., 2011). All of these correction functions 
are equal to or greater than 1. Asc) and f(sc,) are often set to 1 in numerical simulations because 
of the lack of parameterizations and data. Shao and Lu (2000) derived the following simple 
expression for an ideal threshold friction velocity: 


told.) fay Peed 4 A )> (5) 


a a Ss 


where pp (kg/m) is the sand particle density (9,=2560 kg/m°); pa (kg/m?) is the air density 
(pa=1.225 kg/m?); g (m/s?) is the gravitational acceleration (g=9.81 m/s”); and Ay and A, are 
constants (Ay=0.0123 and 4;=1.65x10* kg/s*), which account for the magnitude of the 
interparticle cohesive forces and are empirically determined by a wind tunnel experiment 
(Greeley and Iversen, 1985). Note that we did not measure the sand particle size ds at the 
observation sites. In this study, we used the minimum values of the threshold friction velocity of 
sand particles of size corresponding to a soil particle diameter of about 80 mm. 

2.7.1 Surface roughness effect (shear stress partitioning effect) 


We adapted the surface roughness factor of Raupach et al. (1993) based on drag partitioning to 
include the effects of solid elements (i.e., stones): 


fA) =/(l-o, m, A+A, m, A) 3 (6) 
C 

=£, 7 

B. C, (7) 


where øs is the ratio of the basal to frontal area of the roughness elements; ms is a parameter 
accounting for spatial nonuniformity of surface stress that ranges between O and 1 (1 for 
uniform stress and decreasing as nonuniformity increases); J, is the roughness density of the 
solid elements (Raupach et al., 1993); Js is representing the drag partition of a single roughness 
element; Cr is the drag coefficient for isolated roughness elements; and Cs is the drag 
coefficient for the soil surface. In this study, we used the average roughness density of stones at 
each site for representing the roughness density of the solid elements. For Js, Shao et al. (2015) 
recommended using a value between 100 and 400. In this study, following Tan et al. (2019), we 
set ms=0.5 and £;=100. For obtaining the o;, we used the measured values at each site (i.e., the 
average ratio of breadth to height of the 105 stones). We assumed that o, was constant at a given 
site. 

2.7.2 Effects of soil moisture, salt concentration, and soil crust 

Fécan et al. (1998) parameterized the influence of soil moisture as a function of the clay content 
by fitting an equation to the experimental data for various soils. Shao and Lu (2000) used 
observational data to expand the parameterization of Fécan et al. (1998) for several soil textures 
and obtained the following equation: 
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Cae oo E rae (8) 


l+a(w-w,)’ forw2vw, 


where w (m?/m?) is the volumetric soil moisture; a and b are parameters that depend on soil 
texture; and w, (m?/m°) is the threshold soil moisture value, which depends on the soil texture 
(Kang et al., 2011). Our measured soil moisture during the observational periods was always 
lower than the threshold value. Therefore, we set f(w)=1 in the simulation. We found no salt 
accumulation in the study area, so we did not measure salt concentration and set f(s-)=1. A soil 
crust was sometimes present (Table 1), but we had no method for measuring its hardness. Thus, 
we also set f(scr)=1. 


3 Results 


3.1 Roughness density and roughness length 


Table 2 shows the measurement results for soil texture and the parameters for stones (coverage, 
mean stone height, roughness density, and roughness length) at all sites. The soil texture at all 
sites was loam, except for sandy loam and silt loam at the sites Sub18C and Sub18D, 
respectively. The three highest stone coverages were 57.1% at the Sub18A site (S.MT; Table 1), 
47.2% at the Subl19A site (S.SL), and 39.7% at the Sub19B site (S.SL). The next highest 
coverages were 33.9% at the Sub18C site (N.MT), 28.8% at the Sub18B site (N.SL), 22.4% at 
the Sub14B site (S.SL), and 18.6% at the Main site (N.SL). The stone coverage was only 2.9% 
at the Sub18D site (TD) and 0.0% at the Sub14A (TD), Sub19C (TD), and Sub19D (TD) sites. 
Mean stone height, roughness density, and roughness length generally tended to increase with 
increasing stone coverage. 

We examined the relationship between stone coverage and roughness density based on the 
data of this study and based on the data reported by Marticorena et al. (2006) (Fig. 7); we fitted 
a third-degree polynomial to the data. The roughness density increased gradually when the stone 
coverage was from 0.0% to 30.0%; it increased relatively quickly between coverages of 30.0% 
and 50.0% and apparently increased gradually again when stone coverage exceeded 50.0%. This 
result suggests that it may be possible to estimate the roughness density of stones from the stone 
coverage. At the sites Subl8A, Subl9A, Sub19B, Sub18C, Sub18B, Sub14B, Main, and 
Sub18D, the frequency distribution of the ratio of breadth to height was normal (Fig. 8). 


Table 2 Measurement results for coverage, mean height, roughness density, and roughness length of stones at 
each site 


Site Soil Stone coverage (%) sae cae pa ris Eo 
DEI teure Mean SD Mean SD Mean SD Mean SD (cm) 
Main Loam 18.6 2.95 1.75 0.37 0.39 0.15 0.05 0.001 0.0152 

Subl4A Loam 0.0 — 0.00 — 0.0012 
Sub14B Loam 22.4 3.52 1.81 0.30 0.57 0.23 0.07 0.016 0.0318 
Sub18A Loam 57.1 5.84 1.88 0.36 1.05 0.27 0.34 0.039 0.1569 
Sub18B Loam 28.8 3.26 1.85 0.35 0.48 0.19 0.13 0.017 0.0519 
Subisc Sandy 33.9 4.25 1.74 043 065 037 0.16 0.023 0.0878 
Loam 
Sub18D te 2.9 1.01 1.93 0.41 0.47 0.16 0.01 0.004 0.0011 
Subl9A Loam 47.2 4.13 1.89 0.37 1.01 0.35 0.32 0.031 0.1587 
Sub19B Loam 39.7 4.73 1.84 0.39 0.89 0.30 0.24 0.036 0.1628 
Subl9C Loam 0.0 0.00 — 0.0012 
Subl9D Loam 0.0 — 0.00 — 0.0012 


Note: SD, standard deviation; —, no data available. 
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Fig. 7 Relationship between stone coverage and roughness density from this study (black dots) and the study of 
Marticorena et al. (2006) (open circles) 
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Fig. 8 Frequency distributions of the ratio of breadth to height of stones at the sites Sub18A (a), Sub19A (b), 
Sub19B (c), Sub18C (d), Sub18B (e), Sub14B (f), Main (g), and Sub18D (h) 


Although the coverage, mean height, and roughness density differed greatly from site to site, 
the mean ratio was similar among the sites (Table 2). The roughness correction is related to the 
ratio of breadth to height (i.e., the ratio of the basal to frontal area of the roughness elements 
(os) in Eq. 6); thus, it is possible to use the observed mean ratio in the model simulation because 
the ratio is normally distributed. Previous dust model studies that considered the effect of stones 
on the threshold friction velocity have used a value of 1.00 for this ratio, but in this study, the 
average ratio was 1.84. 


3.2 Observed threshold friction velocity and the effects of roughness density 


Figures 9 and 10 show the time series of friction velocity (green and blue lines), the 
instantaneous threshold friction velocities (red circles), and their temporal means (dashed lines) 
at the sites during 1-5 May 2018 and 27 April-8 May 2019, respectively. Sand saltation was 
observed at all sites in 2019 and at three sites in 2018, but almost no saltation was observed at 
the other four sites in 2018. 

Under the natural field conditions, the instantaneous friction threshold at the sites was not 
constant. However, during both measurement periods, the temporal difference was much smaller 
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than the difference in the mean threshold among the sites. The standard deviation of the 
instantaneous threshold friction velocity around the mean threshold friction velocity value at 
each site was between 0.03 and 0.07 m/s, whereas the standard deviation around the average of 
the mean threshold friction velocity among the sites was 0.15 m/s. The highest mean threshold 
friction velocity (0.71 m/s) was recorded at the site Sub19A, where the roughness density was 
highest (0.320) in 2019. The lowest mean threshold friction velocity (0.23 m/s) was recorded at 
the site Sub19D, where no stones were observed in 2019. 

We examined the relationship between roughness density and mean threshold friction 
velocity at seven sites in 2019 and three sites in 2018 (Fig. 11). Observations were also 
conducted at the sites Subl14A, Sub18A, Sub18B, and Sub18D in 2018, but we could not obtain 
the threshold friction velocity because saltation counts were too few. The mean threshold 
friction velocity tended to increase with increasing roughness density. We divided the dataset 
into three groups according to roughness density: Group 1 consisted of 2019Sub14A, 
2019Sub19C, and 2019Sub19D sites (roughness density=0.000); Group 2 consisted of 
2018Main and 2019Main (roughness density=0.050); and Group 3 consisted of 2018Sub14B 
and 2019Sub14B (roughness density=0.070). The four digits before the site name indicate the 
year of observation. 
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Fig. 9 Friction velocity when the wind direction was between 260° and 370° (i.e., 10°) (blue lines) or between 
10° and 260° (green lines), instantaneous threshold friction velocity (red circles), and mean threshold friction 
velocity (dashed lines) during 1-5 May 2018 at the sites Sub18C (a), Sub14B (b), and Main (c). w+tmean, mean 


threshold friction velocity; SD, standard deviation of the instantaneous thresholds around the mean. 


3.3 Simulated threshold friction velocity and validation against observations 


Figure 12 shows the relationship between the observed threshold friction velocity and simulated 
threshold friction velocity by using the measured roughness density of stones at each site. The 
simulated threshold friction velocity was the same within each of the three groups: Group 3, 
0.42 m/s; Group 2, 0.37 m/s; and Group 1, 0.20 m/s. These values corresponded to roughness 
densities of 0.070, 0.050, and 0.000, respectively. The observed and simulated threshold friction 
velocities were almost the same at the sandy site Sub19D and at the stony sites 2018Sub14B, 
2018Main, Sub18C, Sub19A, and Sub19B (see Table 1), but the simulated threshold friction 
velocities for the sites in groups 1—3 tended to be substantially underestimated. 
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Fig. 10 Friction velocity when the wind direction was between 260° and 370° (i.e., 10°) (blue lines) or between 
10° and 260° (green lines), instantaneous threshold friction velocity (red circles), and mean threshold friction 
velocity (dashed lines) from 27 April to 8 May 2019 at the sites Sub19A (a), Sub19B (b), Sub14B (c), Main (d), 
Sub14A (e), Sub19C (f), and Sub19D (g). 


4 Discussion 


4.1 Soil texture, stone coverage, topography, and dust source 


We found no relationship between the spatial distribution of soil texture and topography, but 
parameters of stones such as coverage and roughness density tended to increase from lowland 
(i.e., TD) to highland (i.e., MT) sites. In addition, those parameters tended to be higher at the 
southern sites (i.e., S.MT and S.SL) than at the northern sites (1.e., N.MT and N.SL), which 
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Fig. 11 Relationship between roughness density and mean threshold friction velocity observed at the bare or 
stony (i.e., no vegetation and no crust) sites (black dots), stony and vegetated sites (open rectangles), vegetated 
site with crust (open triangle), and site with crust only (open rhombus) in spring 2018 and 2019. Symbols show 
the average roughness density and mean observed threshold friction velocity, and bars show their standard 
deviations. For sites Main, Sub14A, and Sub14B, where observations were carried out year-round, the four 
digits before the site name indicate the year of the observations. 
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Fig. 12 Scatter plot between the observed threshold friction velocity (during 1-5 May 2018 and during 27 
April-8 May 2019) and simulated threshold friction velocity using the observed roughness density of stones at 
the bare or stony (i.e., no vegetation and no crust) sites (black dots), stony and vegetated sites (open 
rectangles), vegetated site with crust (open triangle), and site with crust only (open rhombus). The two 
horizontal thin dashed lines show the simulated threshold friction velocity using roughness densities of 0.002 
and 0.020 for barren and sparsely vegetated land cover types, respectively. 


resulted in the threshold friction velocities being low in lowlands and high in highlands. Our 
results agreed with the findings of Callot et al. (2000), who revealed that the size of non- 
erodible elements such as stones was dependent on the distance from the highest topographic 
point. 

It is well known that dust sources are often distributed in topographic lows and lands 
adjoining topographic highs (e.g., Shao, 2008); in this study, such areas are represented by the 
TD sites. Our findings that areas with almost no stone coverage were associated with low 
threshold friction velocities are consistent with this well-known observation about dust sources. 
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4.2 Effects of other factors on threshold friction velocity than stones 


We found that sites with the same roughness density (Figs. 11 and 12) can have different 
observed threshold friction velocities. The roughness density of stones was 0.070 at the sites 
2018Sub14B and 2019Sub14B (Group 3), 0.050 at the sites 2018Main and 2019Main (Group 2), 
and 0.000 at the sites 2019Sub14A, Sub19C, and Sub 19D (Group 1). Increased soil moisture and 
salt accumulation can increase the threshold friction velocity. However, during our observation 
periods, the observed soil moisture was very low (less than 0.05 m?/m°), though it was observed 
at the sites Main, Sub14B, and Sub14A only; further, no salt accumulation was observed. 

In 2019, dead vegetation was observed at the sites Main, Sub14B, and Sub14A, and green 
vegetation was observed at the site Sub14A, whereas no vegetation was observed at these sites 
in 2018 (Fig. 13). Many studies have shown that vegetation reduces wind erosion intensity 
through coverage and shear stress partitioning effects, which decrease soil erodibility, and as a 
result, the sand saltation threshold increases (e.g., Lancaster and Baas, 1998; Li et al., 2008; 
Yan et al., 2011). We can infer that the presence of herbaceous vegetation (live or dead) led to 
higher threshold friction velocities at the sites Main and Sub14B in spring 2019. 


(a) 2018Sub14A 


c) 2019Sub14A Wee dame ten eames [(d) Sub19C 


Fig. 13 Crust conditions at the sites Sub14A (a) and Sub18D (b) in spring 2018, and at the sites Sub14A (c) 
and Sub19C (d) in spring 2019 


In the field, the effects of stones and vegetation on wind erosion cannot be separated from 
those of other factors, such as changes in soil properties or the presence of a crust (Lopez et al., 
1998). We observed crust at the sites Sub14A and Sub18D in 2018 and at the sites Sub14A and 
Sub19C in 2019 (Fig. 13), but the characteristics of the crust differed among the sites. In 2018, 
the observed crust was mostly undisturbed and sand particles were less mobile at the sites 
Sub14A and Sub18D. In 2019, the area of disturbed crust was larger at the site Sub14A, but 
some less-mobile sand particles were still seen, and at the site Sub19C, although a crust was 
observed, many mobile sand particles were present. 

It is very difficult to quantitatively estimate the threshold friction velocity under crusted 
conditions because many factors, including its morphological and compositional properties (i.e., 
the modulus of rupture, thickness, and clay, calcium carbonate, water soluble material, and 
organic matter contents), need to be taken into account (Gillette et al., 1982). However, a well- 
developed crust can protect soil surfaces, thereby reducing wind erosion intensity (Belnap and 
Gillette, 1998). Marticorena et al. (1997a) showed that the threshold friction velocity was much 
higher when crust was undisturbed than when it was disturbed. 
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Crust was observed at the sites 2019Sub14A and Subl19C (Fig. 13). In the plot of the 
simulated threshold friction velocity against the observed threshold friction velocity (Fig. 12), 
the site 2019Sub14A was far to the right of the 1:1 line, indicating that the observed threshold 
friction velocity was much higher than the simulated threshold friction velocity based on the 
observed roughness density of stones alone. Because the site has vegetation as well as crust, we 
cannot isolate the effect of the crust alone, but this result demonstrates that vegetation and crust 
both increase the threshold friction velocity. 

The site Sub19C was also located on the right of the 1:1 line in Figure 12. In addition, the 
observed threshold friction velocity was slightly higher at the site Sub19C than at the site 
Sub19D. The crust may have caused the threshold friction velocity to be higher at the site 
Sub1l9C than at the site Sub19D. However, these sites may be located nearer the 1:1 line 
because of sand particle mobility; many mobile sand particles were observed at the sites 
Sub19C and Sub19D and were relatively fewer at all other sites. 

As discussed in Section 4.1, the topographic depression may be a major dust source because 
stones are almost completely absent there. However, vegetation and crust were distributed 
mostly in the topographic depression (Table 1). Gantsetseg et al. (2017) reported that the 
Normalized Difference Vegetation Index (NDVI) was high in the topographic depression in 
Tsogt-Ovoo during the rainy summer of 2012. Their result supports our finding of dead 
vegetation mainly in the topographic depression in spring because dead vegetation represents 
vegetation that grew in the previous summer or earlier. We could not quantitatively examine the 
effects of vegetation and crust on the threshold friction velocity due to a lack of data, but we 
could determine that an increase in the threshold friction velocity was resulted from their 
effects. In addition, the gaps in the observed threshold friction velocities between sites Sub19D 
and 2019Sub14A, between sites 2018Main and 2019Main, and between sites 2018Sub14B and 
2019Sub14B quantitatively implied the importance of vegetation and crust at these sites. 


4.3 Effects of stones on the threshold friction velocity in the three-dimensional model 
studies 


It is currently difficult to account for the effects of stones on the threshold friction velocity in a 
three-dimensional dust model because no reliable data on the spatial distribution of stones are 
available. However, several studies (e.g., Darmenova et al., 2009; Xi and Sokolik, 2015; 
Foroutan et al., 2017; Lee and Park, 2020) have conducted a three-dimensional dust simulation 
by employing a double shear stress partition scheme. These studies prepared a lookup table that 
contains representative roughness density values of solid elements (i.e., stones) and vegetation 
for each land cover type. The studies used roughness density values of 0.002 (Darmenova et al., 
2009; Foroutan et al., 2017) and 0.020 (Xi and Sokolik, 2015; Lee and Park, 2020) for barren 
and sparsely vegetated land cover types, respectively. In this study, all of the sites were 
considered barren or sparsely vegetated, the threshold friction velocity simulated with 
roughness density of 0.002 was 0.21 m/s, and that simulated with roughness density of 0.020 
was 0.28 m/s (horizontal dashed lines in Fig. 12). Even though all of the sites in this study were 
classified as barren or sparsely vegetated, our results showed that the threshold friction 
velocities simulated using with the observed roughness density varied from 0.20 m/s at the bare 
sites in the topographic depression, to 0.70 m/s at the stony sites on the southern mountain with 
roughness density of 0.000; and these simulated threshold friction velocities were almost the 
same as those observed values at the sites with roughness density of 0.320. 

When we compared the threshold friction velocity estimated by the three-dimensional dust 
models employing a constant roughness density with threshold friction velocities estimated for 
bare or stony sites (black dots in Fig. 12), we found that the threshold friction velocity estimated 
by the three-dimensional dust models employing roughness density of 0.002 was almost same as 
that at a sandy site (Sub19D) in the topographic depression, but the three-dimensional dust 
models employing roughness density of 0.020 overestimated the threshold friction velocity at 
that site. Moreover, the three-dimensional dust models employing both roughness density values 
of 0.002 and 0.020 consistently underestimated the threshold friction velocity at the sites on 
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slopes and mountains (e.g., 2018Main, 2018Sub14B, Sub18C, Sub19B, and Sub19A). Our 
survey result showing that the roughness density varies over a wide range even among sites with 
the same land cover type suggests that reliable spatial distribution data of stones based on field 
measurements are necessary to achieve realistic three-dimensional dust simulation results. 

We expect it to be possible to obtain spatial distribution data of stones over a wide area by 
satellite remote sensing techniques. Some studies have already estimated roughness parameters 
from satellite data. Greeley et al. (1991) found an empirical correlation between the 
aerodynamic roughness length of non-erodible elements and the radar backscatter coefficient 
obtained by Airborne Synthetic Aperture Radar (AIRSAR), and Marticorena et al. (2004) 
obtained a good relation between the aerodynamic roughness length and the Polarization and 
Directionality of the Earth's Reflectance (POLDER)-1 bidirectional reflectance distribution 
function (BRDF). 

Laurent et al. (2005) conducted a three-dimensional dust simulation using a roughness length 
map obtained by the BRDF method of Marticorena et al. (2004), but the mapped data comprise 
roughness values due to both stones and vegetation. For a simulation that discriminates the 
effects of stones and vegetation, such as a simulation using a double shear stress partition 
scheme, satellite observation techniques that can separately quantify stones and vegetation are 
needed. 


5 Conclusions 


We evaluated the effects of stones on the sand saltation threshold for different roughness 
densities of stones by conducting field observations under natural surface conditions during 
sand and dust storms in the northern part of the Gobi Desert, Mongolia. We found that sand 
saltation threshold noticeably increased with the increasing roughness density of stones. 
Specifically, roughness density increased from lowland to highland sites. At the sites to the 
north of Tsogt-Ovoo, where neither crust nor vegetation was observed, the roughness density 
values were 0.000, 0.050, and 0.160 at the TD, N.SL and N.MT sites, respectively. To the south, 
the values were 0.000 and 0.070—0.320 at the TD and S.SL sites, respectively. The stone 
coverage also increased from lowland to highland sites; at the sites to the north, where neither 
crust nor vegetation was observed, the mean threshold friction velocities (the average of the 
instantaneous threshold friction velocities) were 0.23, 0.41, and 0.57 m/s at the TD, N.SL, and 
N.MT sites, respectively. To the south, they were 0.23 and 0.45-0.71 m/s at the TD and N.SL 
sites, respectively. 

We simulated the threshold friction velocity with the model, which includes the Raupach's 
roughness correction that uses measured roughness density values. We found that this roughness 
correction worked very well in the study area. The threshold friction velocity was simulated 
very well at all sites where neither crust nor vegetation was observed, although the simulated 
threshold friction velocity showed a small negative bias. Though our survey was conducted in 
only one small area in the vast Gobi Desert, our results suggest that if roughness density of 
stones over a wide area can be obtained, the roughness correction for the stone effect in the 
model of roughness correction can be applied in the three-dimensional dust models. The 
simulation results of this study also showed that, to reasonably simulate sand saltation in the 
Gobi Desert, it is essential to consider the effects of vegetation and crust as well as the effects of 
stones. 
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